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^ ■ ABSTRACT 

o : 

O ■ We study properties of steady states (states with time-independent 



density operators) of systems of coupled harmonic oscillators. For- 



"^ . mulas are derived showing how adiabatic change of the Hamiltonian 

■*««> . transforms one steady state into another. It is shown that for infl- 

ow, nite systems, sudden change of the Hamiltonian also tends to produce 

5-1 . steady states, after a transition period of oscillations. These naturally 

arising steady states are compared to the maximum-entropy state (the 



thermal state) and are seen not to coincide in general. The approach 
to equilibrium of subsystems consisting of n coupled harmonic oscilla- 
tors has been widely studied, but only in the simple case where n = 1. 
The power of our results is that they can be applied to more complex 
subsustems, where n> I. It is shown that the use of coupled harmonic 
oscillators as heat baths models is fraught with some problems that do 
not appear in the simple n = 1 case. Specifically, the thermal states that 
are though to be achievable through hard-sphere collisions with heat- 
bath particles can generally not be achieved with harmonic coupling 
to the heat-bath particles, except approximately when the coupling is 
weak. 
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I. INTRODUCTION 

Quantum systems with quadratic Hamiltonians have been studied extensively. One reason for this 
is that they are just about the only quantum systems whose time evolution can be found analytically. 
Hence they have provided useful and tractable models in a wide range of fields, including statistical 
mechanics, quantum optics and solid state physics. 

In this paper, we derive some useful formulas for the subset of such systems known as harmonic 
chains. These may be pictured as a chain of coupled harmonic oscillators where the coupling between 
any two oscillators depends only on the distance between them. We derive the most general steady 
states of coupled harmonic oscillators, and then study how adiabatic and sudden changes of the 
Hamiltonian convert one such steady state into another. 

The main motivation for studying harmonic chains has been the pursuit of a dynamical basis 
for equilibrium statistical mechanics. In short, we frequently assume that the density matrix of a 
system is p ex e^^''^^ , i.e. the thermal state corresponding to a temperature T, and it would be 
rather embarrassing if we had no completely solvable heat bath model that explicitly evolved systems 
into such states. An excellent summary of the early developments is given by Huerta and Robertson 
(see [1] and references therein). The main idea is to study the time evolution of the mean values and 
correlations of some small subset of the oscillators, called the system, by taking a partial trace over the 
rest of the oscillators, called the heat bath. The goal is to show that the system will exhibit standard 
thermodynamic features such as Brownian motion and approach to thermal equilibrium. 

Numerous models have been developed that show such features, and we will give only a brief 
overview here. The dynamics of an infinite, one-dimensional chain of oscillators bound only to their 
nearest neighbors was solved in 1914 by Schrodingcr [2]. Klein and Prigogine [3] used these results 
to show that suitable initial conditions for the heat bath led to equipartition of energy in the system. 
Hemmer [4] studied the same model in more detail and found that the heat bath could produce 
Brownian motion of the system. Results for more general quadratic Hamiltonians have since been 
derived by numerous authors [5 — 11]. An up-to-date overview of the models is given by Ford, Lewis 
and O'Connell (see [12] and references therein), who summarize six successful models and show that 
they are all unitarily equivalent. Most of the recent literature on the subject has focused on reduced 
descriptions of a single oscillator in the chain such as a Langevin or Fokker-Planck equation. 

Thus the approach to equilibrium of subsystems consisting of n coupled harmonic oscillators has 
been widely studied, but only in the special case where n = 1. This is rather unfortunate, since as we 
shall see, some interesting complications can arise when n > 1. We will study subsystems consisting 
of an arbitrary number of oscillators, including the seldom mentioned case where the mean values of 
the heat bath oscillators are non-zero. We do this using a model involving cyclic matrices, similar to 
that of [5] . As will be seen below, the great advantage of this model is that subsystems with large 
n can be treated with the same ease as the n = 1 case, since we need not distinguish beetween the 
system and the heat bath a priori. 

This paper is organized as follows: In section II, we establish some notation and derive the most 
general steady state of coupled harmonic oscillators. In section III, we review some useful properties 
of cyclic matrices. In sections IV, V and VI, we study how adiabatic and sudden changes of the 
Hamiltonian transform one steady state into another. Finally, in sections VII and VIII we compare 
these naturally arising steady states with thermal states and discuss some shortcomings of coupled 
harmonic oscillators as heat bath models. A rigorous convergence proof for the sudden change case is 
given in Appendix A. 



II. THE GENERAL STEADY STATE 

In this section, we establish some notation and derive the covariance matrix of the most general 
time-independent state of coupled harmonic oscillators of equal mass. 

As our quantum system, let us take 2N + 1 coupled harmonic oscillators of equal mass, labelled 
—TV, ..., —1, 0, 1, ..., N. Denoting a point in the 2{2N + l)-dimcnsional phase space by 



^1 (1) 

H = — p^p + ^q^Aq, 3 

Im I 

where the time-dependent matrix A is symmetric and positive definite. Throughout this paper, we 
will use units where ni = ujq = \. The number of oscillators can be either finite or infinite. Let us 
define the mean vector /i and the covariance matrix C by 



and the corresponding operators by 



we can write the Hamiltonian as 



fM 



< 



C ^\ I - A^M 

i(qp^ + pq^) (pp^) 



(4) 



(The symmetric ordering is necessary since q and p do not commute.) It is well-known that since fl 
is quadratic, the time evolution of /i and C is identical in classical and quantum statistical mechanics 
and is given by 



{c{t) ^U{t)C{Q)U{tf, 



(5) 
where the time-evolution matrix is defined by 



^t^-\-A{t) OJ^' (6) 

(7(0) = /. 

In the special case where the matrix A is independent of t, this equation can be integrated, yielding ^ 



U{t) = exp 



I\ 


" 




rh 


-A oj 





cos 



A^'H A-^'^sinA^'H 



-Ai/2sinAi/2i cosA^/H 



(8) 



' Here and throughout this paper, the action of a function on a symmetric matrix is defined as 
the corresponding real-valued function acting on its eigenvalues: Since all symmetric matrices A can 
be diagonalized as 

A = RAR^, 

where R is orthogonal and A = diag{di} is diagonal and real, we can extend any mapping / on the 
real line to symmetric matrices by defining 

f{R diag{d^}R^) = R diag{/(d,)}i?^- (7) 

It is easy to see that this definition is consistent with power series expansions whenever the latter 
converge. For example, 

cosAi/2 = yi-i{-A". 

n=0 ^ ' 



One useful feature of ,u and C is that if the phase-space distribution (the Wigner distribution [13]) 
is Gaussian, then they determine the state completely. For all the steady states considered below, this 
will indeed be the case. (Such Gaussian states in quantum optics are often referred to as multimode 
thermal squeezed coherent states, but we will refer to them as Gaussian for brevity.) For a good 
review of Wigner distributions, see [14], [15] and references therein. 

fx and C are also very useful for studying subsystems, as we will do in sections VII and VIII. 
Rather than having to take cumbersome partial traces of the density operator, the means and internal 
correlations of any subsystem can simply be read off as subsets of the entries in /i and C . If the 
state of the total system is Gaussian, so is the state of the subsystem, which means that the latter is 
uniquely determined by these subsets. 

By a steady state of a time-independent Hamiltonian, we mean a state whose density operator 
(or equivalently its Wigner function) is independent of time. The most general steady state is clearly 
a statistical mixture of energy eigenstates, i.e. a density operator of the form 

p = Y,Wn\En){Enl (9) 

where Wn > and ^ Wn = 1. In terms of the Wigner function, this becomes 

W{q,p)=Y,^nW,M,P)- (10) 

Since the energy eigenstates Wn are cumbersome products of Laguerre polynomials, this expression 
does not shed much light on the physical nature of steady states. Fortunately, we can obtain a simple 
and useful result for their covariance matrices: By differentiating Eq. (5) with respect to t using Eq. 
(8), it is readily shown that the mean vector and the covariance matrix will remain constant over time 
if and only if they are of the form 

M =0, 

D \ (11) 



^ ' AD 

where D is a symmetric, positive definite matrix that commutes with A, and A is independent of 
time. We will discuss the case /i 7^ in Section VIII. Note, however, that a covariance matrix of the 
above form will remain constant over time regardless of the value of /i. 

In section IV, we will see how adiabatic change transforms one such steady state into another. 
In section V, we will see that these steady states are a form of attractors, in that quite general states 
tend to converge towards them as i — > 00. In section VI, we will study a one-parameter subset of 
steady states, the thermal states, and compare them with the states arising for adiabatic and sudden 
change. 



III. CYCLIC MATRICES 

In many parts of this paper, we assume that the potential matrix A is cyclic, which means that we 
can write Aij = a\i-j\ and interpret the system as a chain of harmonic oscillators where the coupling 
between any two oscillators depends only on the separation between them. For finite N, we will for 
computational simplicity identify N + i with —N — 1 + i and picture a ring of oscillators rather than 
an array with two ends (thus Aij = Aki ii i — j — k — I mod (2A^ + 1)). Using (7), we can write any 
function of a (cyclic or non-cyclic) matrix A as 



Jy-^Jmn — 



TV 

E 

k=-N 



RmkRnkf{\)- 



(12) 



Cyclic matrices have the great advantage that they all commute. This is because they can all be 
diagonalized by the same matrix i?, an orthogonal version of the discrete Fourier matrix. Physically, 
this means that plane waves form a complete set of solutions. Ford, Kac and Mazur [5] (hereafter 
referred to as FKM) show that if A is symmetric, positive-definite, cyclic and infinite-dimensional, 
then Eq. (12) reduces to 



J\A)mn — 



1 

2^ 



def[\^{e)\ cos{m-n)e, 



(13) 



where the spectral function X^iO) is the function whose Fourier coefficients are row zero of A. The 
spectral function can be interpreted as a dispersion relationship, A being the frequency of a wave with 
wave number 9. Note that f{A) is cyclic as well, i.e. its components depend only on the distance to 
the diagonal. 

A cyclic potential frequently discussed in the literature is the nearest neighbor potential, the case 
where each mass is coupled only to a fixed spring and to its nearest neighbor: 



H 



N 

E 

k=-N 



^Pk + 2«fe + Y (*- 



(14) 



i.e. Akk = 1 + 27^, Ak,k±i = —7^ and all other elements of A vanish. For this special case, the 
spectral function is 



X\0) = 1 + 47^ sin 



2 ^,-2 



(15) 



IV. ADIABATIC CHANGE 

In this section we will study the time evolution of C during an adiabatic change of the Haniiltonian 
of a steady state. 

For generic time-dependent cyclic Hamiltonians, the time-evolution operator U{t) depends on the 
entire history A(t), < t < i, according to Eq. (6). For the case where A changes extremely slowly 
relative to the eigenfrequencies, this simplifies so that U{t) depends only on the values of A at r = 
and T = t, and is independent of A at intermediate times. According to a well-known result, the time 
evolution of a single harmonic oscillator [N = 0) where the frequency lo changes adiabatically is given 
by [16] 

/c.(i)-V2 Wcos0(t) sin0(t)\ /^(0)V2 \ 

^^^>-\ io{tf/^)\-sine{t) cose{t))\ w(0)-i/2j' Ubj 

where 



e{t) = / Uj{T)dT. (17) 

Jo 

We are interested in the case of arbitrary N . Since we are assuming that the matrices A{t) are cyclic 
for all t, they can all be diagonalized by the same time-independent matrix R, 

A{t) = R^ c\i&g{uj,{tf}R, 

where the eigenfrequencies uji also change adiabatically. Hence we can apply (16) to each mode 
separately, which after employing (7) yields 

(A{t)-y^ W cose(i) sine(t)WA(0)V4 \ 

'^^^'-y A{tf'^)\-s{nQ{t) cosQ{t))\ A{Q)-^'^ ) ' ^^^' 



e{t)= / A(T)i/2dT. (19) 



where 



If our system is in a cyclic steady state initially, i.e. if C(0) is cyclic, then D{0) is cyclic and commutes 
with all the above matrices, so substituting (18) into (5) yields 

^w-(T A(t)%t))' (20) 

where 

D{t)^[A{0)A{t)-'Y^'D{0). (21) 

Note that the quantity A(t)D{t)'^ is an adiabatic invariant, i.e. stays constant over time. We see that 
adiabatic change forms an Abelian transformation group on the set of all steady states, and that any 
given steady state can be transformed into exactly those steady states that have the same adiabatic 
invariant. 



V. SUDDEN CHANGE FOR INFINITE SYSTEMS 

In this section we will study the time evolution of C after a sudden change in the Hamiltonian 
of a steady state. We will mainly be interested in the following question: Does C{t) keep oscillating 
forever, or does the state converge towards a new steady state as i ^ oo? 

If the system is in a steady state 

for i < and the potential matrix in the Hamiltonian (3) changes abruptly from Aq to a new constant 
value A ai t — 0, then the time evolution for t > will be given by (5) and (8). Let us temporarily 
choose the simple initial conditions Aq = Dq = I, and return to the general case later. For this choice, 
C(0) = /, so that Eq. (5) gives C{t) = U{t)U{t)^, i.e. 

/ cos2 [A^/H] + A-i sin^ [A^/H] [A-i/2 _ A^/^] sin [A^/H] cos [A^/H] \ 

^^*>-[ [A-i/2 - Ai/2] sin [A^/H] cos [A^/H] cos2 [A^/H] + A sin^ [A^H] J ' ^^^ 

We can separate this C(i) into a time-independent part and an oscillating part with zero time average 
by writing C{t) = Ci + C2(i), where 

_i// + A-i 

^1 -n i+A 



Chit) =i 



(24) 
[l-A-^]cos[2A^/H] [A-i/2 _ Ai/2] sin [2Ai/2i] ^ ^ ' 

[A-i/2 „ ^1/2] sin [2A^/H] [I - A] cos [2^1/2^] 

In Appendix A, we show that for a generic cyclic infinite-dimensional symmetric positive definite 
matrix A, C2{t) ^ as t ^ oo, so that C{t) roughly speaking evolves as follows: 

/ 0\ /Timc-dep.\ f ^^^ \ , .„,. 

l)^[ mess )^[ 1±4 j ^t^°°- (25) 

These results are readily generalized to arbitrary cyclic initial conditions at t = 0: 

X Z\ fD 



^Z Y ) \0 AD' ^«*-°«' (26) 

where 

D^ -[X + A^^Y] . (27) 

If the initial state is the steady state (22), this reduces to 

D = ^[l + AoA-^]Do. (28) 

The fact that all matrices are cyclic is crucial, since it means that they can all be simultaneously 
diagonalized and hence all commute. 

Note that whereas adiabatic change formed a transformation group on the set of steady states, 
sudden change (and subsequently waiting until the oscillations have died down) does not. The steady 
state reached through a sequence of two sudden changes, one after the other, can in general not be 
achieved with a single sudden change. This is obvious from the fact that no more than half of the 
energy can be lost in single sudden change. 

A second difference from the adiabatic case is that sudden changes in general have no inverse, i.e. 
cannot be undone. If the potential is suddenly changed from Aq to A and then, after the oscillations 
have died down, back to ^Oi the net result is 



D^-[2I + A^^A + AoA-^]D. (29) 



1 

4 
This steady state has a greater energy than the initial one for any A ^ Aq. 
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VI. SUDDEN CHANGE FOR LARGE BUT FINITE SYSTEMS 

The result C2{t) ^ as t ^ oo does not hold for finite N. Indeed, for finite N, the components 
of 6*2 (t) must return to their initial values an infinite number of times, the typical time between such 
recurrences being the Poincare recurrence time. How can this be reconciled with the results of the 
previous section? An excellent discussion of these matters is given by Mazur and MontroU [6] for some 
special cases. Here we will limit ourselves to an informal heuristic discussion, whose aim is to give a 
qualitative feeling for the non-cyclic case and the N ^ oo limit. 

Applying Eq. (12) to the components of the various terms in C2{t), the function /(A^) always 
contains a trigonometric factor that oscillates increasingly rapidly with respect to A^ as t -^ oo. 
Suppose for a moment that the product RmkRnk varies smoothly with k and that the spectrum \\ 
is not too degenerate. Then for very large t, all phase coherence between different terms in the sum 
is lost and we can for all practical purposes replace the arguments of the trigonometric functions by 
random variables 4>i that are uniformly distributed on the interval [0,27r]. Using {cos (pi) = and 
{cos (f)i con cf)k) ~ ^Sik, this means that if we let t be a random variable uniformly distributed on the 
interval [— T, T] for some large enough T, then the expectation value and variance of a typical term 
like 

X{t) = Acos\2A^/h] (30) 

is approximately given by 

N N 

{X{t),nn) ^ {Y^ i?„fei?„fcA^ COS (2Afet)) = Y^ i?™fci?„feA^(cOS</>fc) =0, (31) 



k=-N k=~N 



N N -, W 



(A(i)f„„) = Y Y RmkRnkRmlRnl><i>H {cos 4>k cos (j)l) ^ - ^ {RmkRnk>'k) ■ (32) 

k=-Nl=-N k=-N 

Since the typical element of a generic N x N orthogonal matrix R is of order N^^/"^, the sum above 
will be of order N x A^~^, which yields a standard deviation of order N^^^^. Thus for large but finite 
systems, we expect the elements of C2(t) to evolve as follows: 

(i) During an initial transition period whose duration is of the order of the dynamical time scale 

uJq , they decay from their initial values of order unity to values of order N~^''^ . 
(ii) After that, they oscillate around zero with an oscillation amplitude of order N^^^"^. 
(iii) Since their time evolution is almost periodic, they must return to values of order unity an 
infinite number of times. This happens approximately once every Poincare recurrence time. 
However, as shown by Mazur and MontroU, the Poincare time scale is generally enormous 
compared to the dynamical time scale, since it tends to grow exponentially with N for systems 
of this type. 



VII. IMPLICATIONS FOR HEAT BATH MODELS 

In this section we will discuss thermal states, and compare them to the stationary states arising 
from adiabatic and sudden change. It will be seen that the latter two are in general not thermal. 

Thermal states 

The energy of a steady state is 

Given a Hamiltonian and a fixed energy, the thermal state is the state that has the maximum entropy 
consistent with this energy. The Wigner function of a thermal state is Gaussian. FKM show that a 
thermal state of temperature T is the steady state given by 

= 

_ 1 /A-i/2coth/3Ai/2 \ , (34) 

^ 2 (^ Al/2coth/3Al/2 j 

where we have defined /3 = -^^ . In the classical limit /3 <C 1 , this reduces to 
whereas T = yields the ground state 

a.i(V ;.)■ <-' 

Heat bath models 

According to the standard interpretation, a system will be in a thermal state if it has been exposed 
to an ideal heat bath long enough to reach thermal equilibrium. (After this the heat bath can be 
removed if the system is kept isolated, since the thermal state is a steady state.) Numerous attempts 
have been made to model ideal heat baths by infinite systems of coupled harmonic oscillators. We will 
now use the results we have derived for adiabatic and sudden change to illustrate some shortcomings 
of such models, and clarify what they can and can not do. 

Typically, the time evolution of some small subset of oscillators, called the system, is studied by 
taking a partial trace over {i.e. ignoring) the rest of the oscillators, called the heat bath. A goal of all 
such pursuits has been to devise interaction Hamiltonians A{t) such that the state of the subsystem 
becomes thermal. 

Our use of /i and C has the advantage that we need not distinguish between system and heat 
bath a priori: if we decide to view some subset of oscillators as the system, then their means and 
internal correlations are simply given as subsets of the entries in fi and C. 

We can summarize our previous results for steady states as follows: 

(^0^^^) D in adiabatic case, 

i [/ + AqA"^] D in sudden case, (37) 

i^-i/2coth/3^^/^ in thermal case. 

From this it readily follows that we can transform an arbitrary steady state given by Aq and Dq into 
a thermal state with an arbitrary temperature T by adiabatically changing Aq into 

2 




A = 



icoth-i(2Aj/2z?o 



(38) 



However, we are not interested in such "fine-tuned" models, where A{t) is chosen to depend on the 
initial data Dq. Rather, we are looking for a fairly robust heat bath model, where one single interaction 
potential can thermalize fairly general initial states. 



Initially uncoupled oscillators 

As a simple example, let us choose the initial state to be the thermal state at some temperature 
To for a Hamiltonian given by Aq = I. Since this corresponds to uncoupled oscillators, this state is 
simply a direct product of single-oscillator states: /i = and C = ^ coth(3ol, where Pq = j^- Now 
let a cyclic interaction potential A be switched on, cither adiabatically or suddenly. After a while 
(after the switch-on is complete in the adiabatic case or after the oscillations have died down in the 
sudden case), we can compare the resulting steady states with the thermal state: 

/^-l/2 Q \ 

^" ^ M n 41/2 ) coth/3o after adiabatic change. 



C= <^ 



I+A- 







Cs = ^ I 2 ^ J coth /3o as t ^ 00 after sudden change, (39) 

_ 1 /A-i/2coth/3yli/2 

^*~2l Ai/2coth/3yli/2 



for a thermal state. 



We see that these three covariance matrices are in general all different. In particular, the final states 
in the adiabatic and sudden cases are not thermal. Although there is no reason to expect the total 
correlation matrices to be thermal, one might hope that certain submatrices of Ca or Cs should be 
thermal, in the sense that they equal the corresponding submatrices of Ct — at least if Ct is allowed 
to refer to a different temperature, i.e. allowing T 7^ Tq. In other words, one might hope that certain 
subsystems would approach the state that they would have if the total system were thermal t. It 
should be fairly obvious that this is not the case for generic A. The nearest-neighbor coupling is a 
simple counterexample. For this simple case, we can evaluate the momentum correlations explicitly 
in the limit 7^1, kT 3> IIluq: 



(PmPn) 



j^Tc;^^ ^_^^^_^^^^ after adiabatic change, 

l^cothA, ^ [(1 + 7')<5™,„ - ^(<5™,„-i + <5™,„+i)] after sudden change, (40) 

TT^Sjn n in thermal case. 



(The adiabatic result has been derived using Eq. (13) and 7^1, the sudden result is exact and the 
thermal result is only valid for kT 3> fiwo-) We see that the qualitative behavior is quite different in all 
three cases: In the adiabatic case, weak long-range correlations exist. In the sudden case, correlation 
exist only between nearest neighbors. In the high-temperature thermal case, there are no correlations 
at all. 

For |?Ti — n| <C 7, the integral (13) for position correlations is completely dominated by the region 
around = 7, and we obtain 



ln(27r7) 
27r7 



after adiabatic change, 
(qmqn) ~ { \^m..n + ^ after suddcu change, (41) 

-r^-a in thermal case. 

47/3 

Note that for the adiabatic and thermal cases, this is independent of m and n , which physically means 
that oscillators with separation |?7i — n| ^ 7 move as one rigid unit, having their positions almost 
perfectly correlated. 



^ This is perhaps the only reasonable definition of what to mean by a subsystem becoming thermal- 
ized. Recall that we have only defined thermality of a system relative to some Hamiltonian. Thus if a 
subsystem in some sense becomes thermal, the question is: Relative to which Hamiltonian? It would 
be naive to expect the answer to be the Hamiltonian describing only the forces within the subsystem, 
since when a subsystem is coupled to its surrounding, the Hamiltonian that governs it by definition 
involves parts outside of the subsystem. 
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The simple 2x2 submatrices corresponding to a single oscillator are all diagonal, just as they 
are for the thermal state of a single oscillator. Since they contain only two non-zero elements, they 
are indeed identical to the thermal state of a single oscillator with some temperature T and coupling 
uieff 7^ ti>, or identical to the 2x2 submatrices of Ct with some temperature T and coupling jeff 7^ 7- 
In this limited sense, both adiabatic and sudden change can be said to thermalize n = 1 subsystems. 
For subsystems with n oscillators, Ca, Cg and Ct each have n(n+ 1) different non-zero elements in the 
submatrix, so we see that the above-mentioned match by adjusting the two parameters T and Jeff 
was just a fortuitous coincidence that worked when n — I. 

Another popular choice of initial conditions in the literature has been a thermal state with some 
coupling Aq 7^ /. Also this case can be readily treated with the above formalism, and generally fails 
to yield thermal states for n > 1 subsystems. 

Does this mean that cyclic quadratic Hamiltonians are totally incapable of acting as heat baths 
for n > 1, or is there some simple limit in which subsystems can be thermalized? One frequently 
discussed limit is that of weak coupling. 

The weak coupling limit 

In the case of infinitesimally weak coupling, say A = I + eF, Equations (39) reduce to 



/ — |F \ coth/3o 



C= <^ 



I+^F) 2 

I — jF \ coth 00 



I+^F) 2 



after adiabatic change, 

as t ^ oo after sudden change, (42) 



where 



^ [1 + «] 2^ \ ££tM for a thermal state, 

1 + [I - a\ ^t J ^ 



2/3 



sinh 2/3 



(43) 



Thus we see that for infinitesimally weak coupling, the sudden case coincides with the adiabatic case, 
but that the resulting state is generally still not thermal. It should come as no surprise that the 
sudden case coincides with the adiabatic case in this limit, since a sudden infinitesimal change is in a 
sense adiabatic. We see that there are only two cases when the resulting state is thermal: 

The first is if F is a multiple of /, which gives T ^ Tq and corresponds to the trivial case of the 
oscillators remaining uncoupled. 

The second is when a = 0, which implies T = Tq ~ 0. This simply reflects the fact that adiabatic 
change changes the ground state into the new ground state. Note that this is a pure quantum 
phenomenon, since a^l7^0as?i^0. 

The sense in which harmonic heat baths work 

Although the above results for the weak-coupling limit showed that subsystems did not become 
thermal, they did indeed become almost thermal, the difference being terms of order e. This is not 
very impressive at all, in view of the fact that the initial state was also almost thermal. However, 
a general feature of infinite harmonic chains is that the final state after sudden change remains the 
same even if the initial states of a finite number of oscillators is changed. This follows directly from 
the fact that all components of U{t) approach zero as i -^ oo. This means that a sudden change will 
result in the same almost thermal state even if some subsystem starts out in a completely different, 
quite non-thermal state. It is in this sense that the many harmonic heat bath models in the literature 
can transform a subsystem into an approximately thermal state. 

Unfortunately, the interesting case where the system has internal couplings stronger than those 
of the heat bath cannot be treated with the formulas derived in this paper, since the matrix A would 
not be cyclic. It would appear, however, that the strength of the coupling between the system and 
the heat bath must be carefully balanced for such a heat bath model to work. If the coupling is so 
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strong that it is comparable to the couplings within the system, then the final state of the system will 
probably deviate considerably from the thermal state. If, on the other hand, the coupling is too weak, 
then the relaxation time, the time scale over which the system approaches the thermal state, will be 
unreasonably long. 

VIII. THE INDEPENDENCE OF DISPLACEMENT AND DISPERSION 

In this section, we will discuss how the previous results for sudden change are modified in the case 
of non-zero means, i.e. when /z 7^ 0. It will be seen that the fact that /i and C evolve independently 
of each other makes systems with /i 7^ totally unable to produce thermal states. 

We can generalize (3) to arbitrary quadratic Hamiltonians by writing 

H = z'^ Bz = tr Bzz^ , (44) 

where B is any symmetric matrix. This covers both cases where the different particles have different 
masses and cases where the potential is not positive definite. 

(H) = tr B(zz^) = tr B^ifi^ + tr BC = ^i + E^, (45) 

say, where we will call Ei = tr B/j./^^ the displacement energy and and E2 = tr BC the dispersion 
energy. (To be precise, the second equal sign is valid only if the two off-diagonal n x n submatrices of 
B are symmetric and thus identical — otherwise the non-commutativity of q and p gives rise to an 
uninteresting extra additive constant.) We see from (5) that fj, and C evolve completely independently 
of each other. By formally setting either /x or C equal to zero, (5) thus shows that the displacement 
energy and the dispersion energy are separately conserved, independently of one another. Hence 
displacement energy can never be used to increase dispersion and entropy. This is in stark contrast 
to an approach to the thermal state, where the resulting equilibrium entropy for subsystems is the 
maximum allowed by the total available energy. This difference between the thermal state and the 
actual state as t ^ 00 is illustrated in Figure 1 . The Gaussian Wigner distribution in the phase-space 
of a single oscillator is shown by the contour at which it has dropped to half its maximum value, a 
circle. Recall that the entropy of the oscillator depends on the area of the circle, whereas the energy 
depends on the average distance to the origin. Thus the circle symbolizing the equilibrium state has 
the maximum area consistent with the available energy. The actual state has the same energy but 
much less entropy. 

Thus oscillator systems which are able to evolve a state with non-zero means into something 
resembling the thermal state cannot have purely quadratic Hamiltonians, but must contain some 
non-linear couplings. 
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IX. CONCLUDING REMARKS 

We have studied steady states of coupled harmonic oscillators, and how adiabatic and sudden 
changes in a cyclic potential transform one steady state into another. Our main conclusions are as 
follows: 

(i) Adiabatic change can transform a given steady state into any other steady state, subject to 
the constraint that the adiabatic invariant AD^ is conserved. 

(ii) After a sudden change, the covariance matrix undergoes an oscillatory phase, after which 
it generally converges to a new steady state. This asymptotic behavior depends crucially 
on how degenerate the spectrum of A is. Generic sudden changes are not reversible, and a 
sequence of sudden changes, where the final and initial potentials A are the same, leads to a 
net increase in energy. 

(iii) The states resulting from adiabatic and sudden changes are in general not thermal states, 
not even if we restrict our attention to subsystems. We conclude that harmonic heat bath 
models can lead to approximately thermal subsystems only under the following conditions: 

1) The coupling between the system and the heat bath is weak. 

2) /i = 0, i.e. the expectation values of the displacements and momenta of the oscillators 

vanish. 
We summarise the state of affairs for heat bath models as follows: It is frequently assumed that 
the density matrix of a system is p ex e^^'*^^, i.e. the thermal state corresponding to a temperature 
T. Unfortunately, no completely solvable heat bath model has ever been found that explicitly evolves 
multiparticlc {n > 1) systems into such states. The problem appears to arise when the coupling 
between the system particles and individual heat bath particles is so strong that it is comparable to 
the couplings within the system. Thus approximately thermal states can be created only when the 
when the system particles are very weakly coupled to the bath particles. If this is the case, the system 
particles must be coupled to very many heat path particles to avoid the approach to equilibrium taking 
an unreasonably long time. Only in the limit when a typical system particle becomes infinitesimally 
coupled to infinitely many bath particles (such as in the limit of the FKM model [5] or the 10 model 
[12], both dealing with the n = 1 case) is the limiting state exactly thermal. Thus we may expect 
nature to be full of approximately thermal states, but devoid of exactly thermal states. 
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APPENDIX A 

In this appendix, we will discuss the circumstances under which C2(t) — > as i ^ oo. 

Applying Eq. (12) to the components of the various terms in C2(i), the function /(A^) always 
contains a trigonometric factor that oscillates increasingly rapidly with respect to A^ as i ^ oo. In 
Section VI, we gave a qualitative discussion of how this causes C2(t) to approach zero for general 
symmetric positive definite A with non-degenerate spectra, as well as what happens when N is large 
but finite. Here we will give a rigorous treatment for the special case when A is infinite-dimensional 
and cyclic. 

The behavior of C2{t) depends crucially on how degenerate the spectrum of A is. Since all cyclic 
Hamiltonians are invariant under parity, they all have the degeneracy 

uj{-e) = uj{e), (46) 

where we have defined the frequency lo = |A|. If this is the only degeneracy, i.e. if w{9) is invertible 
for > 0, and if in addition lo{9) is differentiablc, then we can change variables in Eq. (13) and obtain 

/(A)„„ = 2/ f{io^)g{uj)cos[{m-n)e{cu)]duj. (47) 



Here we have introduced the spectral density function g{w), defined so that g{u})cLj is the fraction of 
the frequencies that lie in the interval [u',u' + cLu]. Note that 

ZTT dlO 

if 9{ijj) is single- valued and invertible, but that g{uij) is well-defined as a distribution for any spectrum 
uj{6). For the case of nearest-neighbor coupling, 

{ -—, = if 1< Iwl < 1 + 472, 

1 otherwise. 

We wish to examine the various terms in C2 [t) . Let us define 

O-n = J\A)m,m+n (5U) 

and choose the function to be 

f{A) = cos{A^/h). (51) 

Eq. (47) now yields 

Gnit) — 2 / cos{ujt)g{uj) COS [n6{uj)]dLj. (52) 



By extending the functions g and 9 to the entire real line by g{—Lo) = g{^), 0{^^) ~ ^('^)i ^^d 
g{u}) ~ unless LOmin ^ l^^l ^ ^maxi wc scc that this can be written as a Fourier transform: 

an{t) = Un{t), (53) 

where 

Un{oj) = g{uj)cos[n9{uj)] . (54) 

According to Riemann-Lebesgue's Lemma, 

Un{t) -^ as |i| ^ 00 (55) 
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if Un{uj) is an integrable function. Since the latter is trivially satisfied for all n (indeed J_ Un{uj)duj ~ 
27r(5on), this implies that ^ 

cos(Ai/2i) ^0 as t ^ oo. (56) 

In an almost identical fashion it is readily shown that sin(A^'^f) ^ as t ^ oo. By multiplying these 
results by various constant matrices, it follows that all terms in Ciit") approach zero. In summary, we 
have shown that all components of C2(i) ^ as t ^ oo for any symmetric, positive definite, cyclic, 
infinite-dimensional A if its spectral function A^(0) is invertible and differentiable for > 0. 
It is straightforward to generalize this result in a number of ways: 

* If A^(0) is differentiable but not invertible because it turns around a finite number of times, 
then the same result is obtained by integrating over each monotonic segment separately. 

* If A is merely non-negative definite, then A^^/^ and A~^ are undefined. However, the result 
still holds for the momentum-momentum part of the covariancc matrix. 

Roughly speaking, the key result C2(i) -^ as i — > oo fails if the spectrum is too degenerate. 
An easy way to appreciate this is to study the diagonal terms of the matrix cos(A^/^i), which simply 
equal 

ao(i)=.9(t). (57) 

For uncoupled oscillators, we have A oi I and the spectrum is totally degenerate, i.e. the spectral 
density gii^^ = 6(lu — luq) + S{ll! + luq) and ao(t) oc cosLUot for all t. By Riemann-Lebesgue's Lemma, 
doit) ~^ as t -^ oo if the spectral density is an integrable function, not if it is merely a tempered 
distribution such as 6. Furthermore, from partial integration of the Fourier integral, we know that if 
the spectral density is a k times differentiable function, then an{t) approaches zero faster than t^*^. 



T The convergence is not necessarily uniform, i.e. although cos(^^'^t)m„ -^ for any fixed m 
and n, we do not always have sup„ „ cos{A^''^t)mn ^ 0. A more careful analysis shows that uniform 
convergence is obtained if the spectral function A is analytic on the entire interval [— tt, tt] and in 
addition is nonlinear. 
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